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1. Straight Forward 

Goals: We strive for Matlab fluency on 4 planes; To speak the language 
free of syntatic error, to confidently translate English simulation and/or design 
questions, to identify the proper numerical method for the job, and to represent 
the solution in the most visually striking way possible. 

Means and Evaluation: We believe these goals are best attained by im- 
mersion, i.e., by regular exposure to fluent Matlab speakers in a number of 
distinct engineering settings. Though distinct, these settings will all reduce to 
problems of network simulation or design and will introduce you to a number 
of important numerical methods. The fluent speakers are the Rice Learning 
Assistants that will be conducting the numerous small sections of CAAM 211. 
In order to bring you to fluency we have constructed weekly two part Matlab 
excursions. 

The first part is a short Owlspace Quiz over the science, math and program- 
ming constructs of the week. These will be "multiple choice," not timed, and in 
"working" each quiz you may consult written and electronic media but you may 
not discuss the quiz with any living person. They are due early each Friday and 
will be reviewed in your Friday lab. The quiz and lab are stuctured in order to 
facilitate your solution of the weekly project. 

The second part is a full MATLAB project where you will be asked to trans- 
late English instructions into a MATLAB program that produces visual output. 
For two Assignments you will pledge that you received no aid in their comple- 
tion. For the other assignments I encourage you to compose and troubleshoot 
with one another, but require you to write, document, festoon and submit your 
own work. By document I mean that you write your own English preamble 
(header) and that you write your own English comments throughout the code. 
By festoon I mean that you must label axes and title plots in an original and 
creative fashion. The main header must conform to this standard: 

First Name Last Name, CAAM 210, Fall 2010, HW N 
file name 

program description 
program usage instruction 
program usage example 

Here is one example, we shall see many more. 
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While you are encouraged to work with other students currently registered 
in CAAM 210, you are not allowed to consult with students who have taken 
previous offerings of CAAM 210 for assistance with homework assignments (ex- 
cept the 211 instructors and designeted college labbies). Use of external tutors 
must be cleared through your 210 instructor. The appearance of code obviously 
originating from previous semesters is considered to be a serious honor code 
violation. 

You will publish each assignment to a single pdf file that you will subse- 
quently upload to our Owlspace page. I will now explain how to publish. 

• Create the directory CAAM210. 

• Save this mfile (as type trigplot.m) to your CAAM210 directory. 

• Fire up Matlab, and set your current directory to CAAM210, by typing 
"cd CAAM210" at the Matlab prompt. 

• Type "edit trigplot" 

• From the Editor Window click on "File" and select "Publish Configura- 
tion for trigplot.m" and follow the menu down and right to "Edit Publish 
Configuration for trigplot.m" 

• This will bring up an "Edit M-File Configurations" window. Look for "Out- 
put settings" at about midwindow. Click and hold on the "html" as "Output 
file format" then drag down and select "pdf or "doc" if you are using an 
older version of MATLAB. 

• In the "Output folder" row please type "CAAM210" (without quotes) 

• Finally click on "Publish" at the bottom right and check for the appearance 
of trigplot.pdf in your CAAM210 directory. 

The pdf will contain a copy of your code, a list of things written to the 
command window, and each of the figures your code generated. If you are not 
happy with your pdf please close the pdf reader before trying again. 

Grades: Your grades will be determined solely on the basis of the 14 quizzes, 
the 14 assignments and your lab participation. 

211: Your score for CAAM 211 will be (Q + P)/2 where Q is your average quiz 
score (out of 100) and P = 100a/12 where a is the number of the 13 Friday lab 
sessions that you attended. Absences may be excused and late quizzes may be 
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accepted only with written confirmation by your doctor, academic advisor or 
college master. 

210: Your score for CAAM 210 is 

( e i + e2 + X>) 

where pj is your score for project j (each out of 100 except the two pledged 
projects which are each out of 200) and is your score on the kih extra credit 
essay (each of which is out of 50). Late projects and essays may be accepted only 
with written confirmation by your doctor, academic advisor or college master. 

Resources: Make use of your instructor's and/or teaching assistant's office 
hours. If these hours are inconvenient please ask to schedule an appointment. 
The graduate teaching assistant offers an optional help/review session every 
Thursday 7-9pm. The Friday lab is staffed by a trained Rice Learning Assistant. 
It is not optional. In the lab you will work in small groups of 3 to 5 students. 
The lab will not be equipped with computers. Please bring your laptop. The 
packet you are now reading is the textbook for this course. The definitive 
Matlab resource resides online at Mathworks. We will make constant use 
of this. I strongly recommend that you take advantage of the deep student 
discount and purchase a copy of Matlab from Mathworks. It costs $99 at the 
Mathworks Store. 
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2. Loops, Matrices and Fractals 

To iterate is to apply the same rule over and over again. Though simple to 
say and easy to code this can produce some amazing patterns, e.g., eye balls, 
palm trees and galaxies. 

To get started consider the rabbit rule 'multiply by 6'. If we start with a 
population of 10 and iterate we achieve 10, 60, 360, 2160, 12960, ... No, I am 
not amazed by this pattern. If we instead apply the cannibal rule 'multiply by 
1/2' we arrive at 10, 5, 2.5, 1.25, 0.625, ... again, nothing new here. 

To get more interesting patterns we move into higher dimensions. More pre- 
cisely, instead of transforming one number into another via scalar multiplication 
we wish to transform a pair of numbers, say x(l) and x(2) : into a second pair, 
say y(l) and y(2), via matrix multiplication as depicted below. 






Figure 2.1. A pair of vectors in the plane. 
More precisely, we write 



y(l) = A(l,l)*a;(l) + A(l,2)*ar(2) 
y{2) = A(2, 1) * x(l) + A(2, 2) * x{2) 




A4(l,l) A(l,2)\ 
VA(2,1) ,4(2,2); 
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For example, the matrix 



A = 



( 



1 
-1 



) 



rotates every vector clockwise by 90 degrees. Do you see it? Try it by hand 
for a few vectors and then coax Matlab to do it by trying variations on the 
following sequence of commands: 

>> A = [0 1; -1 0] ; 

» x = [1; 1]; 

>> y = A*x; 

» plot ( [0 x(l) ] , [0 x(2) ] ) ; 

>> hold on 

» plot ( [0 y (1) ] , [0 y (2) ] , ' r' ) ; 

>> y = A*y; 

» plot ( [0 y (1) ] , [0 y (2) ] , ' r— ' ) ; 

>> y = A*y; 

» plot ( [0 y (1) ] , [0 y(2) ] , ' r- . ' ) ; 

>> axis equal 

I agree that repeated application of A will just take us in circles. You yearn 

for something different. 

Visualization in Matlab 

Although MATLAB has a friendly interactive face its power lies in its ability 
to chew on user generated programs. Yes, a program is just a list of legal 
Matlab commands like our little example above. We will be building and 
changing programs of greater and greater complexity and so it is much wiser to 
create and save such programs than to reenter them interactively at the MATLAB 
command line. 

One creates and saves programs through the use of a text editor. MATLAB also 
has a built-in editor that may be invoked it by typing edit from the MATLAB 
prompt, >>. No matter how you create any Matlab program you should save 
it with a '.m' extension, e.g., caam210.m is OK. For this example, typing omg 
at the prompt (assuming that Matlab's current directory contains a program 
omg.m) will execute the program omg.m. 

Here is a program, Circle Deform I, that plots the unit circle and its deforma- 
tion under a particular matrix transformation. This program waits till the end 
to show what it has done. To watch it work on the fly check out Circle Deform II 
and Circle Deform III. These programs demonstrate the use of a for loop. They 
also sport an informative header. They also enjoy proper indenting and are 



5 



sprinkled with many comments. They each produce a plot with an informa- 
tive title. All of these elements, as well as the if clause are in evidence in 
our Circle Deform IV. You will use a for loop and an if clause in your first 
project. 

I recommend that you save the Circle Deforms to your directory, run them, 
change the A matrix, run them, change A, run them until your excitement 
ebbs. These runs are static incarnations of Matlab's dynamic eigshow demo. 
I encourage you to run the demo, hit the help button, and ponder the significance 
of eigenvalues and eigenvectors. 

An eigenvector of a matrix A is a vector that gets stretched by A but not 
rotated, and the eigenvalue is the amount of stretch. The geometry is clear, but 
is there any value in pursuing these notions? Eigen is German for self. Ask 
Google to Google itself by searching for eigenvalue and google. Or better yet, 
check out the $25,000,000,000 Eigenvector 

Project: Fractal Fern 

We will be growing ferns by random repetition of a few simple matrix/vector 
transformations. We have seen, on the lecture page, how to multiply a 2-by-2 
matrix and a 2-by-l vector and we have interpreted the "action" via stretching 
and rotation. A third natural geometric transformation is to shift or trans- 
late. This is done algebraically by simple vector addition. In particular, if 
x=[x(l);x(2)] and y=[y(l);y(2)] are two 2-by-l vectors then their sum z=x+y 
is also 2-by-l and its elements are z(l)=x(l)+y(l) and z(2)=x(2)+y(2). 

The old-school fern below was produced by this code. You see that at each step 
it generates a random number and applies one of two possible transformations, 
the first with probability 0.3 and the second with probability 0.7. 

Your task is to dissect and document this existing code and then to write 
and run code (called fern.m) that produces a significantly more interesting fern. 
In particular, your new program should successively generate a random number 
and should 

apply z = [0 0;0 0.16]*z with p = 0.01 

apply z = [0.85 0.04; -0.04 0.85]*z + [0; 1.6] with p = 0.85 

apply z = [0.2 -0.26; 0.23 0.22]*z + [0; 1.6] with p = 0.07 

apply z = [-0.15 0.28; 0.26 0.24]*z + [0; 0.44] with p = 0.07 
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Figure 2.2. A fern assembled from two afflne transformations. 

You should accomplish this with a single if clause, although, given that 
your logic now forks four ways, you will want to make use of some subordinate 
else if statements. 

Once you are happy with your work please publish a single pdf via following 
the instructions in the Straight Forward. 
Your work will be graded as follows: 

1 figure: 10 pts. Use MATLAB title command to place your name on plo 

1 M-file: 10 pts for header 

10 pts for further comments in code 

5 pts for indentation in for loop 

5 pts for indentation in if clause 
10 pts for use of elseif 

Your header must conform to the standard in the Straight Forward. For more 
headers and examples of inline documentation return to our circle deformers. 
Please recall that each individual is solely responsible for commenting and 
headering their code and for labeling their own plots. 

For those who wish to dig deeper, please check out the text, Fractals Ev- 
erywhere by M.F. Barnsley, or hit the wiki page on Iterated function system. 



3. Solving one Real Equation via Bisection 

Matlab provides many means with which to solve equations. The diary 
demonstrates the use of the polynomial root finder, roots, and the symbolic 
solver, solve. In response to your question, "Hey man, how do these solvers 
solve?" we shall devote the rest of the semester to writing Matlab solvers of 
increasing complexity and hence applicability. 

For real scalar equations with real solutions there is a simple divide-and- 
conquer algorithm that is fun to visualize and easy to code. Let us take, for 
example, the function 

f(x) = x 5 + x 2 - 2 

and attempt to solve f(x) = 0. We start by judicious snooping and note that 
/(0)/(3/2) < and so f(0) and /(3/2) have opposite sign and so the Interme- 
diate Value Theorem permits us to conclude that a solution lies in between 
and 3/2. To bisect now means to check the sign of /(0)/(3/4) and to restrict 
attention to the interval [0,3/4] if negative or to [3/4,3/2] if positive or to stop 
if /(3/4) is small enough. This process is repeated until the value of / is indeed 
small enough. Here are the main steps. We suppose we have a Matlab function 
that evaluates f(x) for a given x. We also suppose the user has provided us 
with an interval of interest, [a, b] and a tolerance t. Our job is to find an x in 
[a, 6] for which \f(x)\ < t. 

1. If f(a)f(b) > inform the user that a solution may not exist in [a, 6] and 
exit. 

2. Set x = (a + b)/2. While \f(x)\ > t do 

3. If f(a)f(x) < set b = x else set a = x. Go to step 2. 

In order to code this you shall build a function that (i) accepts the arguments, 
a, b, and tol, then (ii) executes steps 1-3 above and finally (hi) returns the 
answer. 

In order to get you off on the right foot let me provide the following skeleton 
(it performs only 1 bisection) and a diary of its use 

>> x=biskel (0, .5,1) 

Sorry, but I can not be sure your fun changes sign on [a,b] 
x = NaN 

>> biskel (0,1.5,1) 
x = 0.7500 
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I have left it for you to code steps 2 and 3. There are a number of ways to do 
it. I suggest a while loop that encloses an if clause. 



Project: Bar Cooling 

We consider a bar of length L. We suppose its temperature is zero at one end 
and, at the other, is proportional to the outflux of heat there. The rate at which 
the bar dissipates heat is the square of the least, strictly positive solution, x, of 

sin(xL) + xcos(xL) = (3.1) 

Our task is to determine how this root depends on the length of the bar. In 

particular, we will reproduce the figure below. 

Steve Cox - Assignment 2 




1.5 



2 2.5 3 

Length, L 



3.5 



Figure 3.1. The decay rate is a decreasing function of bar length. 

We will accomplish this by composing 3 simple functions within a single m-file. 
The m-file (called cooldrive.m) will look like 

function cooldrive % this is the driver 
% set a, t and a range of L and b and find x 
% and plot x"2 against L 
return 



function x = cool (a, b, t , L) 

% code bisection 

return 



% for a given a, b, t, and L find x 



function val = coolfun (x, L) 
val = sin (x*L) + x*cos (x*L) 
return 



% for a given x and L evaluate "cool" 
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In cooldrive I would set a = 0.1 and 6 = 3 and t = 0.01 and compute the root 
x at each L = 1 : .1 : 4 and produces a figure like that above. I recommend 
accumulating the L and x values and issuing one plot command at the end of 
the loop. For example, 

for cnt = 1:31 

L (cnt) = 1 + (cnt-1) /10; 
x(cnt) = cool (a, b, t , L (cnt ) ) ; 

end 

plot (L, x . "2) 

The trouble with this loop however is that (3.1) may have more than one root 
between a = 0.1 and b = 3. I recommend that you move one or both of these as 
you change L. 

Your work will be graded as follows: 

The M-file: cooldrive. m is a function that calls cool, 
a function that calls coolfun. 

15 pts for detailed headers of cooldrive, cool, and coolfun 
featuring 'usage' and 'examples' and definitions 
of all arguments and outputs (as in biskel.m) 

10 pts for further comments in code 
5 pts for indentation 

10 pts for correct code. 

Plot: 10 pts for plot of x"2 versus L labeled as above but with your name. 
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4. Solving one Real Equation via Newton's Method 

Although bisection is simple and easy to code there exist far better methods 
(faster and applicable in more situations) for solving nonlinear equations. The 
best methods stem from Newton's application of what is today called 

FTE : The Fundamental Trick of Engineering: Pretend the world 
is linear. 

More precisely, suppose / is a real function of the real variable x, e.g., f(x) = 
x 3 — 8, and we wish to solve f(x) = 0. We begin with the hope that the solution 
is near some point xq. For points x near xq we recall that Taylor's Theorem 
permits the expression 

f{x) = f(x ) + f(x )(x - x ) + f(x )(x - x ) 2 /2 + f'"(x )(x - xof/S\ + • • • 

where • • • stands for higher order terms, meaning terms like (x — xq) p where 
p > 3. Now Newton's Method consists in applying the FTE to the above as 
a means of providing a better guess. Namely, let X\ be the place where the 
linear part of / near xq vanishes. That is, solve 

= f(x ) + f'(x )(xi - Xq) 

for x\. That is, set 

X\ = Xq - f(x )/f'(x ) 

If / was indeed linear we would be done, i.e., / would vanish at x\. In the 
nonlinear case X\ is not a solution but is rather a (hopefully) better guess than 
what we started with. So you now wonder, "Can this new guess be improved?" 
The answer is yes. It does not require any new ideas, just stubborness. Namely, 
apply the FTE (repeatedly) until you arrive at a guess you deem "close enough" . 
Anything you do repeatedly cries out to be looped. The fundamental Newton 
Step to be taken each time through the loop is 

X J = X J-l ~ f( X j-l)/f( X 3-l) 

One should exit the loop when is sufficiently small. 

For a slideshow graphical interpretation of each Newton Step hit this Newton Demo 
site. 

Here is a diary of Newton's Method applied (successfully) to f(x) = x 3 — 8 
starting from xq = 1. 

Here is a diary of Newton's Method applied (unsuccessfully) to f(x) = x 3 — 
2x + 2 starting from xq = 1. 
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Project: Newton vs. Bisection 

Get to know this Deluxe Bisection routine. Write a similar Deluxe Newton 
routine, i.e., one whose first real line is 

function [x, iter] = denewt (x, t , L) 

and that also solves coolfun (x, L) =0 for x when given L. 
Write a driver function (solve race) that sets 

L = 1, a = 0.1, 6 = 3, x0=(a + 6)/2, 

calls deluxe bisection and denewt (with starting guess xO), over the range of 

tolerances t = 10 - - 7 , j = 1,2,..., 8 and produces the figure below. You will 

need to ask MATLAB for help with legend and semilogx. 

Bisection vs. Newton 

30 1 1 1 i I 




Tolerance 

Figure 4.1. The lesser your tolerance the more you should prefer Newton. 

Your code should reside in a SINGLE mfile called solverace . m and it 
should look like 

% solverace header 
function solverace 

code calls debis and denewt and plots their iters per tol 

return 

% denewt header 

function [x, iter] = denewt (x, tol , L) 

code Newton' s method, calls coolfun and coolfundx 

return 

% debis header 

function [x, iter] = debis ( a, b, t ol , L) 
code Bisection, calls coolfun 
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return 

% coolfun header 

function val = coolfun (x,L) 

code — evaluate coolfun 

return 

% coolfundx header 

function val = coolfundx (x, L) 

code — evaluate the derivative, with respect to x, of coolfun 

return 

Your work will be graded as follows: 

solverace, 6 pts for header 



6 pts 
2 pts 
8 pts 



for further comments 
for indentation 
for correct code. 



in code 



denewt , 



6 pts 
6 pts 
2 pts 
10 pt 



s 



for header 

for further comments 
for indentation 
for correct code. 



in code 



The plot, 



4 points 
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5. Solving one Complex Equation via Newton's Method 

We now show that Newton's method extends easily to one complex equation 
in one complex unknown. For example, let us see whether Newton can confirm 
that the roots of 

z 2 - 2z + 2 

are in fact 

1 + i and 1 — i. 

The Newton rule remains 

z = z- f(z)/f(z) 

and so I hope it is obvious to you that a real starting guess will produce only 
real iterates and will not get us back to our complex roots. Please open my 
meandering diary and then run a few examples for yourself. Now, with a nonreal 
guess, our diary indeed has a happy ending. Now, if all we wanted were the 
roots, we would simply type roots ( [1 -2 2] ) at the MATLAB prompt. Our 
interest however, is not the destination but the journey. Along the way we will 
gain an appreciation for the limitation of Newton's method and we will grow 
our pallete of Matlab art tools. 

In particular, we will dissect Isaac's attack on a polynomial by first assigning 
each root a color. Then, depending on where Isaac takes us we paint accordingly. 
For example, if Isaac takes us to the first root we paint it red, if he takes us to 
the other root we paint it blue, if Isaac loses his way we paint the point black. 
As he has a much better chance of losing his way in the complex plane we will 
produce a number of gorgeous Newton wastelands. 

You will write code that dissects a user specfied quartic. To get started we 
dissect the particular 

(z 2 -l)O 2 + 0.16) (5.1) 
Although the main Newton step is of course 

z = z - (z 2 - l)(z 2 + 0.16)/(4z 3 - I.682) (5.2) 

it is unclear which, if any, root this may lead to. It all depends on where we 
start. We shall see now that (5.2) beautifully partitions the complex plane into 
5 distinct regions, four so-called Newton Basins and one remaining Newton 
Wasteland. The Newton Basin associated with the root z = 1 is defined to 
be all those starting points for which (5.2) eventually produces z = 1. The 



14 



other basins are denned accordingly. The Newton Wasteland corresponds to all 
those starting points for which (5.2) fails to converge to a root. For more details 
and examples please browse the Student Gallery and hit this Newton Basin site. 
With respect to our example, we choose the following color scheme 

Basin of z=l is yellow 
Basin of z=— 1 is red 
Basin of z=0.4i is green 
Basin of z=-0.4i is blue 
Wasteland is black 

With this scheme this pointilist driver produces the picture below 

for x = . 15 : . 0025 : . 55, 
for y = -.15: .0025: .15, 

color = newt (x+i*y, le-3, 20) ; 

plot (x, y, color) 

hold on 

end 
end 




Figure 5.1. A sampling of the Newton basins of (5.1). 

Let us reverse engineer what newt must be up to. Well, it takes a seed x + iy 
and a tolerance (le-3) and a maxiter (20) and returns a color ('y.', 'r.', 'b.', 
'g.', or 'k.') depending on where (5.2) took the seed. The driver then paints 
that seed accordingly. The range of x and y specify the window of the complex 
plane to be painted, while the increment (0.0025) specifies the resolution. Here 
is the code that does the job. 
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Though this is indeed straightforward to understand, it does not exploit Mat- 
LAB's ability to work directly on vectors and matrices. We shall see that the 
nested for loops can be contracted to a single (very fast) line. 

x = . 15 : . 0025 : . 55; 
y = -.15: .0025: .15; 
[X,Y] = meshgrid (x, y ) ; 
Z = X+i*Y; 
for k=l rmaxiter, 

Z = Z - (Z. *2-l) .* (Z. ~2+0.16) ./ (4*Z. ~3 - 1.68*Z); 

end 

Here Z denotes the full grid of complex Newton seed's, and, thanks to the 
dot operator (recall quad) we may operate on all seed's simultaneously. 

Now, upon completing this for loop, Z will be a matrix whose elements are 
either close to one of the four roots, or not. To find out which ones are close to 
1 I recommend the find command 

[il,jl] = f ind (abs (Z-l) <0 . 1) ; 

Now you simply place yellow at (x(jl),y(il)). In order to see what is hap- 
pening please consult this 9 seed diary. If we paint the other roots accordingly, 
and leave the wasteland blank, we arrive at the finer (a markersize of 1) portrait 
below. 

z 4 -0.84z 2 -0.16 




Figure 5.2. Newton basins of (5.1). 

Our assignment this week is to produce such portraits for arbitrary quartics. 
Quartics are encoded by 5 complex scalars, for example, the quartic 

2z A - 4z 3 + (2 - i)z 2 - iz + 10 
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is simply 

q = [2-4 2-i -i 10] 

In order to evaluate a polynomial at a given z one just types 

f = polyval (q, z ) 
Now, the derivative of the polynomial above is simply 
dq = [8 -12 4-2i -i] 

And to evaluate it at a given z one just types 

df = polyval (dq, z ) 

You may use polyval for this assignment but I ask that you write your own 
polyder. You should give it a new name and test it against polyder. 

The last innovation for this assignment is the translation of the q into a string 
that will be used to title your portrait. Typing help strfun will lead to hours 
of string fun - after which you'll know all you need for this week. But just to be 
sure, the line 

qlab = strcat (qlab, ' +' , num2str (q (k) ) , ' z " ' , num2str (5-k) ) 

ought to be helpful. Of course indiscriminate use of this will lead to titles like 

lz A + Oz 3 + -0.84z 2 + Oz + -O.I62 

rather than the slick one that adorns our finer portrait. I know I can count on 
you to suppress leading Is, kill terms with leading 0s, and never write +- or z . 
You might use strrep to replace these unsightly strings. 

Project: Newton Basins 
Write a function called qnewt that takes as arguments 

q, a vector of 5 complex coefficients of a quartic 
xt, a vector of 3 x grid values, xlo, xinc and xhi 
yt, a vector of 3 y grid values, ylo, yinc and yhi 
maxiter, the max number of iterations 

and then paints the portion of the complex plane (where xlo < x < xhi and 
ylo < y < yhi) with points colored by root. Your qnewt must follow the 
meshgrid and find procedure sketched above. In using find you will need 
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to ask which elemenents of Z are close to the kth root. The kth root is simply 
r (k) if you have executed r=roots (q) . 

Your qnewt must also make use of your homegrown polyder function. Your 
qnewt must automatically label its portrait with the user defined polynomial, 
as in Figure 5.2, achieved via 

qnewt ([1 -0.84 -0.16], [.45 .0001 .55], [-.05 .0001 .05], 20) 

Drive this code with a driver, called qdrive, that sets the grid parameters as 
in the example above, and calls qnewt on the above quartic and its 3 neighbors 

[1 -0.84 -0.1 -0.16] 
[1 -0.1 -0.84 -0.16] 
[1 -O.li -0.84 -0.16] 

Your work will be graded as follows 

10 pts for headers 

8 pts for further comments in code 

6 pts for indentation 

10 pts for correct myownpolyder 

10 pts for correct qlab (no Is or 0s or +-) 

8 pts for 4 plots, each titled with its associated quartic. 
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6. Boolean Gene Networks 

Newton's Method can be viewed as a Discrete Dynamical System in that at 
prescribed times (corresponding to iterations) we move from point to point in 
the complex plane under the action of a specific rule, z = F(z) where F(z) = 
z — f(z)/f'(z) is a rational function built from the complex polynomial /. In 
fact, what we have done is to consider 

z, F(z), F(F(z)), F(F(F(z))), ... 

either until a new iterate produces no change, because f(z) « 0, or until we 
reach a preordained maximum number of iterations. We shall, in this chapter, 
widen our scope beyond the complex rational F of the last chapter. Now z will 
be the state of a gene network and F will implement the inexorable logic of 
gene regulation. 

Our first goal will be to model and study networks of genes along the lines of 
the work of Wuensche. We shall assume that a gene is either on (1) or off (0) 
and that a gene network is merely a Boolean (after George Boole) Network, i.e., 
a List of nodes, with a wiring list and a logic rule for each node. 

It will be convenient to assume that each gene is regulated by exactly three 
of its neighbors. As such there are then 8 possible inputs at each gene and 
so 2 8 = 256 possible logic rules at each gene. For example, the 6 node net 
associated with Fig 12 in Wuensche is 




Figure 6.1. A network of six genes where each gene is regulated by three 
others. 

This was produced via the network visualization command biograph, included 
in the Bioinformatics Toolbox, from information specified in the wire matrix. 
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wire = [4 2 3; 5 3 2; 3 6 1; 5 4 6; 6 12; 3 5 6]; 
n = size (wire, 1 ) ; 
a = zeros (n + 1 ) ; 
for i=l : n 

a (wire ( i , : ) , i ) = 1 ; 

ids{i} = num2str(i); 

end 

ids {n+1} = 'Gene Net ' ; 
g = biograph (a, ids ) ; 
selfcon = f ind (diag (a) ==1 ) ; 
for i=l : length ( self con) 

g . nodes ( self con ( i )). Shape = 'diamond'; 

end 

view (g) 

As biograph shuns self-regulation we have adopted the convention that such 
genes will be depicted as diamonds. We next express the action at each gene 
via a table that dictates the local logic, 

Table 6.1 

111 110 101 100 Oil 010 001 000 rule 

1 1 1 1 1 1 231 

01000000 64 
1 1 5 

1 1 1 1 108 

1 1 1 1 1 61 

1 1 1 1 1 62 

Several of these logic functions are fairly complex while several are pretty simple. 
For example, gene 2, inhibits itself and requires activation from both genes 5 and 
3. In other words, if s(i) and ns(i) are the respective current and subsequent 
(next) states of gene i then 

ns(2) = (s(5) AND s(3)) AND (NOT s(2)). 

Similarly, the logic at gene 6 follows 

ns(6) = (s(3) XOR s(5)) OR (s(6) AND NOT (s(3) OR s(5))). 

Let us compute a number of exact cases. For example, if 
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s = (10 10 1) 

then upon consulting wire (1, : ) and s we find that 001 is the pattern pre- 
sented to gene 1. Consulting the associated column in row 1 of Table 6.1 we 
conclude that ns(l) = 1. 

Next, on consulting wire (2, : ) and s we find that 010 is the pattern pre- 
sented to gene 2. Consulting the associated column in row 2 of Table 6.1 we 
conclude that ns(2) = 0. 

Next, on consulting wire (3, : ) and s we find that 111 is the pattern pre- 
sented to gene 3. Consulting the associated column in row 3 of Table 6.1 we 
conclude that ns(3) = 0. 

This is indeed very mechanical and so is best turned over to MATLAB. In 
this case we arrive at 

ns = (1 1). 

As it happens ns is also the next state of itself and so is called a point attractor 
of the gene net described by Table 6.1. 

We now come to the question of how best to diagram the transitions between 
each of the 2 n states of a gene net with n genes. We begin by numbering the 
states from to 2 n — 1 by equating them with there decimal equivalents, e.g., 

(1 1 1) =2~5+2~3+2~0=41 and 
(10 1)= 2~5 + 2~0 = 33. 

We then construct a 2"-by-2 n State Transitition Matrix, STM, of zeros, with a 
single 1 in each row, where STM (i, j ) =1 if state i-1 transitions to state j-1. 
(We subtract 1 because MATLAB row and column numbers begin with 1, not 
0). Once STM is built we may hand it to biograph to produce the associated 
State Transition diagram. 
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18 State Transition Diagram 



41 49 55 




Figure 6.2. The State Transition Diagram for the net specified in Table 6.1. 

We see indeed that state 41 — > 33 — > 33 as computed above. We also see two 
(more fascinating) basins of attraction associated with the 5-state attractor 

50 -> 35 -> 37 -> 5 -> 15 50 

and the 8-state attractor 

1 -> 43 -> 52 



38 7^ 11 -»• 54^ 39 -> 1. 



In this week's project we will investigate the sensitivity of these attractors to 
mutations in the logic rules. 



Write a 



Project: Gene Networks I 



function STM = genestm (wire , rule ) 



that accepts an n-by-3 wire matrix and an n-by-1 rule vector and returns the 
associated State Transition Matrix. 

After determining n, the number of genes, I would translate the rule vector 
into an n-by-8 rule matrix (without the top and right borders) like that in Table 
6.1. I would do this one row at a time via a subfunction of the form 
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function b = d2b(r,C) 

that converts decimal numbers to C-bit binary. I ask that you not use the 
builtin functions bin2dec and dec2bin, for they are too slick for us. Hint: 

b (C-f loor (log2 (r) ) ) =1, reset r, repeat. 

In order to build STM I would proceed like 

for i=l : 2 " n 

s = d2b ( i-1 , n) ; 

ns = next state (using s, wire and rulemat) 

j = b2d(ns); your binary to decimal converter 

STM ( i , j+1) = 1; 

end 

Finally, include your genestm function under a function called genestmdriver 
that takes no arguments, but sets wire and rule (to the values used in our exam- 
ple above), draws the gene net using biograph as above, calls your genestm 
and plots your State Transition Diagram. Next, flip one bit in your rule vector, 
e.g., rule (2) =192 or rule (4) =44, and repeat all of the above steps. 

Your work will be graded as follows: 

10 pts for headers for genestmdriver, genestm, d2b and b2d 
8 pts for further comments in code 
4 pts for indentation 

8 pts for correct rule matrix (displayed to command window) 
8 pts for correct computation of next state 
4 pts for correct STM 

2 pts for one gene net plot 

6 pts for two State Transition Diagrams 
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7. Probabilistic Boolean Gene Networks 

Our experimental understanding of the logic of gene regulation is rarely (if 
ever) as clean as that presented in the last section. It is more often the case that 
regulation follows multiple possible courses. Probabilistic Boolean Networks 
(PBN) were invented to deal with this scenario. In this case, each gene i has 
R(i) rules 

r(i,l), r(i,2),...,r(i,i*(i)) 
and associated probabilities 

p(i, l),p(i,2), . . . ,p(i,R(i)). 

For example, lets consider the fully coupled 3-gene net 




Figure 7.1. A network of 3 genes with 3-way regulation, 
with wire, rules and probabilities 

Table 7.1 

gene wire rule probability 

1 1 2 3 238, 230 3/5, 2/5 

2 1 2 3 64 1 

3 1 2 3 232, 128 1/2, 1/2 

As the number of rule choices varies with each gene we need a new data structure 
for rule. We recommend the cell, e.g., 

rule = { [238 230] 
[182] 

[232 128] } ; 
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With n genes we now have 



N = ]\R{i) 



i=i 



possible (old school, nonprobabilistic) gene networks. In our example case we 
depict the rule tree that leads to these N = 4 gene nets. 



rule 1 





Rule Tree 



rule 2 (1) 



rule 2 (2) 




Figure 7.2 The rule tree for the PBN specified in Table 7.1. The branch prob- 
abilities were displayed via biograph (ptree, ids , ' showweight s ' , 'on') 

We begin with rule 1 at the top (root) of the tree and draw two branches 
(with probabilities) associated with the 2 rule choices. At the next generation 
we choose the rule 2 with probability 1. As rule 3 has two possiblities we again 
branch at each node and arrive at 4 possible nonprobabilistic gene nets. Let us 
develop this in somewhat great detail. In particular, let us build ptree, the 
adjacency matrix that we hand to biograph that then illustrates our rule tree. 

There is first the question of size. How many nodes are there in the rule 
tree? The numerical answer is Ntree = 1 + 2 + 2 + 4. How can you 
build this number automatically from the R vector, where R ( i ) = number of 
elements of rule{i}. Once Ntree is found we can initialize via ptree = 
zeros (Ntree+1 ) and label the first node via ids{l} = 'rule 1'. 

We next specify the first row of ptree by "connecting" the top node to the 
next two nodes with the proper weights 

ptree(l,2) = prob{l}(l) 

ptree (1,3) = prob{l}(2) 
or, in one line, via ptree ( 1 , 2 : 1+R ( 1 ) ) = prob{l}. Next we set their 
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associated labels 

ids{2} = ' rule 2 (1) ' 

ids{3} = ' rule 2 (2) ' . 
We next connect these nodes to the next level of the tree via 

ptree(2,4) = prob{2} ptree(3,5) = prob{2} 
and set their associated labels via 

ids{4} = 'rule 3 (1)' ids{5} = 'rule 3 (2)'. 

Finally we connect to the last level via 

ptree(4,6) = prob{3}(l) ptree(4,7) = prob{3}(2) 

ptree(5,8) = prob{3}(l) ptree(5,9) = prob{3}(2), 

and set their labels via 

ids{6} = 'Net 1' ids{7} = 'Net 2', 

ids{8} = 'Net 3' ids{9} = 'Net 4'. 
This results in the ptree matrix 
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Your challenge is to automate the construction of this matrix, and to assemble 
the associated probabilities and rule vectors of the final N nets. 

The probability, P(i), of arriving at the zth net is simply the product of the 
branch probabilities that connect the zth leaf to the root. For example, the net 
probabilities for the PBN specified in Table 7.1. are 



P(l) 


= prob{l}(l)prob{2}(l)prob{3}(l) = 


3/10 


P(2) 


= prob{l}(l)prob{2}(l)prob{3}(2) = 


3/10 


P(3) 


= prob{l}(2)prob{2}(l)prob{3}(l) = 


2/10 


P(4) 


= prob{l}(2)prob{2}(l)prob{3}(2) = 


2/10. 


For the zth net we may compute the State Transition Mat 


rix, ST Mi, as in the 
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previous chapter. The full Probabilistic State Transition Matrix is then simply 

N 

PSTM = P(i)STMi. (7.1) 

i=l 

3 State Transition Diagram 




Figure 7.3. The State Transition Diagram for the PBN specified in Table 
7.1. 

We observe that states and 7 are point attractors while 6 — > 4 — > 2 — ► 6 
and 6 — > 5 — > 6 are pseudo-attractors. 

Project: Probabilistic Boolean Gene Networks 

Write a 

function [pnet, rnet] = ruletree (rule, prob) 

that accepts rule and prob cells, plots the associated rule tree (see Figure 7.4) 
and returns a vector, pnet, of N net probabilities and a vector, rnet, of N 
rule indicators. 

Append this function to a function called pbndriver that takes no argument 
but sets wire, rule and prob to 

wire = [4 2 3; 5 3 2; 3 6 1;5 4 6; 6 1 2; 3 5 6]; 
rule = {[231 230], 64, [5 7], 108, 61, [62 60]}; 
prob = {[1/2 1/2], 1, [1/2 1/2], 1, 1, [1/2 1/2]}; 

calls ruletree and proceeds to build the Probabilistic State Transisition Ma- 
trix of (7.1) by using last week's code to compute the State Transition Matrices 
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of the N oldschool nets encoded in rnet. Finally, use biograph to depict the 
final State Transition Diagram as in Figure 7.5. 

Regarding the format of rnet I recommend using powers of 10. In this case, 
the output of ruletree should be 

pnet = 0.1250 0.1250 0.1250 0.1250 0.1250 0.1250 0.1250 0.1250 
rnet = 111111 111112 112111 112112 211111 211112 212111 212112 

in which case we see that each net is equally likely, that the first net uses the 
first rule at every gene, that the second net uses the first rule at gene 1 through 
5 and the second rule at gene 6, etc. 



rule 5 (1) 



Net 1 



rule 6 (1) 
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Net 7 





Net 8 



Figure 7.4. Rule Tree for the PBN specified above. 
Your work will be graded as follows: 

10 pts for headers for pbndriver and ruletree 

8 pts for further comments in code 

4 pts for indentation 

8 pts for correct pnet computation 

8 pts for correct rnet computation 
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4 pts for correct assembly of PSTM 

4 pts for ruletree plot 

4 pts for State Transition Diagram 



PTslSt bte Transition Diaqrfc 




Figure 7.5. The Probabilistic State Transition Diagram. 
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8. Stochastic Simulation of Reaction Networks 

Our Boolean approach to gene regulation invoked the logic rule of each operon 
without questioning the underlying biophysical mechanism. The key players are 
the gene (stretch of DNA) and its reader (RNA Polymerase). Governing genes 
regulate the transcription of their subjects by either promoting or prohibiting 
the binding of RNAP to the subject's stretch of DNA. 

This suggests that we might benefit from a tool that simulates the interactions 
that occur in a network of reacting chemical species. In the next two chapter 
we will consider two standard means for modeling and simulating such systems. 
The first is associated with the name of Gillespie and the second with the names 
of Michaelis and Menten. 

Stochastic Chemical Kinetics 

This is a beautifully simple procedure for capturing the dynamics (kinetics) 
of a randomly interacting (stochastic) bag of chemicals (us?). 

We begin with N reacting chemical species, Si, S2, . . . , Sjy, and their initial 
quantities 

We suppose that these species interact via M distinct reactions 

Ri, R2, Rm 

and that these reactions occur with individual propensities. 

Ci,C 2 , ...,c M 

For example, if R\ is 

X\ + X2 — > X3 

then we suppose that the average probability that a particular S1S2 pair react 
within time dt is c\dt. If, at time t, there are X\ molecules of Si and X2 molecules 
of S2 then there are X1X2 possible pairs and hence the probabilty of reaction 
R\ occuring in the window (t, t + dt) is X\X2C\dt. (Under the assumption that 
molecules of Sj are hard spheres of diameter dj and mass rrij, and that the 
reaction takes place in region of volume, V, and temperature, 0, only when the 
associated kinetic energy of the colliding spheres exceeds the activation energy, 
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u\, it follows that the propensity 

-K (di + d 2 \ 2 Ske(m 1 + m 2 ) , * ma \\ 

where k is Boltzmann's constant). 

For more general reactions we will write h m for the number of distinct R m 
molecular reactant combinations available in the state (Xi, X n ). We then 
arrive at the reaction probability a m = h m c m . We shall have occasion to call on 

ao = cli + CL2 H h 

Given this list of reactions and their propensities one naturally asks 
Which reaction is likely to happen next? 

and 

When is it likely to occurr? 

Daniel Gillespie answered both questions at once by calculating 

P{T, m)dT, 

the probability that, given the state (Xi, X n ) at time t, the next reaction 
will be reaction m and it will occur in the interval (t + T, t + T + dT). 

A careful reading permits us to write P(T, m)dT as the product of two more 
elementary terms, namely, the probability that no reaction occurs in the window 
(t,t + T), and the probability that reaction m occurs within time dT of T. We 
have already seen the latter, and so, denoting the former by Pq(T), we find that 

P{T,m)dT = P (T)a m dT 

Regarding Po(T), as the probabilty that no reaction will occur in (t,t + dT) is 
1 — aodT it follows that 

P (T + dT) = P (T)(l-a dT) 

or 

(P (T + dT) - P (T))/dT = -a P (T) 
and which, in the limit of small dT, states 

iftT) = -a P (T) 

and so 

P (T) = exp(-a T) 
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It now comes down to drawing or "generating" a pair (T, m) from the set of 
random pairs whose probablity distribution is 

P(T,m) = a m exp(-a Q T) 

Gillespie argues that this is equivalent to drawing two random numbers, T\ and 
r2, from the [0,l]-uniform distribution and then solving 



for m. We have now assembled all of the ingredients of Gillespie's original 
algorithm: 

Gather propensities, reactions, initial molecular 
counts, and maxiter from user /driver . 
Set t = and iter = 0; 

while iter < maxiter 



Plot each X_j at time t 

Calculate each a_j 

Generate r_l and r_2 

Solve (8.1) and (8.2) for T and m 

Set t = t + T 

Adjust each X_j according to R_m 
Set iter = iter + 1 



end 

The gathering of reactions, like the gathering of wire and rule, may take some 
thought. For small networks we can meet the problem head on. We reproduce 
two examples from Gillespie. The first is 

X^ + Y ^2Y with propensity c\ ^ ^ 

2Y — ► Z with propensity c 2 

Here the subscript on indicates that its level remains constant (either be- 
cause it corresponds to a constant feed or it is in such abundance that the first 
reaction will hardly diminish it) throughout the simulation. Here is the code 
that produced the figure below. 





for T and solving 



ai + a 2 + ... + a m -i < r 2 ao < a\ + a 2 + ... + a 



m 
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Gillespie (29): x+y->2y->z, c x = 5, c 2 = 0.005 



13 
O 

o 

O 2000 



time 



Figure 8.1. The evolution of Y following (8.3). 

Though poorly documented you can nonetheless see each step of the algo- 
rithm. Our second example is 

with propensity c\ 

with propensity c 2 (8-4) 
with propensity C3 



^00 + ^1 
Y x + Y 2 

Y 2 



2Yi 
Z 



Roughly speaking, the Y\ species procreates, is consumed by the Y 2 species, 
which itself dies off at a fixed rate. Such models are known as predator-prey 
and, one can imagine, that the two populations often keep one another in check. 

For example, here is the code that generated the lovely figure below. 

Lotka via Gillespie 
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Figure 8.2. The evolution of Y\ and Y 2 following (8.4). 

Of course for larger reaction networks we, I mean you, must proceed more 
systematically. On the road to a more general robust method we wish to address 

1. The user is unsure of how to choose maxiter and would prefer to provide 
the final time, tfin. 
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2. The user may enter the reaction list as a cell and code the update of h in a 
subfunction and so preclude a very lengthy if clause. 

3. A single run of Gillespie is rarely sufficient - for it delivers but one possible 
time evolution. The user would prefer that we made nr runs and collected 
and presented statistics across runs. 

We may incorporate the first point by exchanging Gillespie's for for a while. 

To address the second point, I imagine a cell where each entry encodes a full 
reaction, say, e.g., by adopting the following convention 

rtab { k } ( 1 : 2 : end) = indicies of reaction species 
rtab { k } ( 2 : 2 : end) = actions for species above 

for the kth reaction. In the Lotka example, this would read 

rtab ={[11] % Yl -> Yl + 1 

[1 -12 1] % Yl -> Yl - 1 and Y2 -> Y2 + 1 
[2 -1] } % Y2 -> Y2 - 1 

and the update function would look like 

function h = update (y) 
h ( 1 ) = y(l); 
h(2) = y(l) *y(2) ; 
h(3) = y (2) ; 

The use of update is relatively straightforward, e.g., 

a = c. * h; 

With this a in hand you may now use cumsum and find to find the next 
reaction index (in just a couple of lines - and with NO if clauses). With this 
index in hand you may visit the proper element of rtab. 

To address the third point above, we should suppress plotting on the fly and 
instead return the time and desired solution vectors to a driver that will collect 
and plot statistics. In particular, if we lay our Gillespie runs along the rows of 
a big matrix and take means (averages) down the columns, e.g., 

for j=l : nr 

[t,x] = mygill (tf in, rtab, xO, c) ; 
X(j, :) = x; 

end 

avg = mean (X, 1 ) ; 
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we "ought" to be on the right track. The trouble is that each of our Gillespie 
runs are loyal to particular time vectors. One way around this is to legislate a 
uniform time vector 

tvec = : tine : tfin 

and bring the individual x vectors into X only after interpolating them with 
interpl. 

Project - Stochastic Operon Simulation 

We will follow McAdams and Arkin in their application of Gillespie's method 
to a simple model of gene expression. In particular we model the binding of RNA 
polymerase, R, to the gene promoter segment, Pr. The complex, RPr, may 
either decompose into its constituents or it may be transcribed and translated 
in a process that creates 10 new molecules of protein, P, in addition to the 
free constituents, R and Pr. Each protein molecule undergoes degradation and 
reversible dimerization (protein-protein binding). We depict these 6 reactions 
via 

R + Pr — > RPr RNAP promoter binding and opening 
RPr — > R + Pr RNAP promoter closing and unbinding 
RPr — ► 10P + R + Pr protein production 

(8-5) 

P + P — > D dimerization 

D — ► P + P dedimerization 
P — > protein degradation 

and collect the reactants into the vector 

x=[R Pr RPr P D] 

and follow their evolution from the initial levels 

xO = [10 1 0] 

and 6 propensities 

c= [2 1 4 2 0.5 0.05]. 

Working from the outside in, you will write a function called gilldriver 
that sets 

tfin, the duration of the simulation 
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tine, the increment between interpolated time points 

nr, the number of Gillespie runs 

rtab, a cell array that encodes the M reactions 

xO, the initial levels of x 

c, the propensities 

and call your subfunction my gill nr times and plots (as below) the mean 
(blue) and plus/minus one standard deviation (red) of the dimer solutions. Your 
my gill subfunction should behave like 

function [t, dimer] = mygill (tf in, rtab, x, k) 

and should call on its associated update function and should not require any 
if clauses. 




time 



Figure 8.3. Unchecked dimer growth. 

In the absense of inhibition the dimer count is roughly increaing. To complete 
this assignment I ask you to augment your gilldriver to incorporate repression 
of Pr by D via the two additional reactions 

Pr + D — > / D binds Pr and so R can not read 
j p r _^ jj pj unbinds Pr 



(8.6) 



I would start from 1 = and assume propensities of 0.01 and 0.005 and produce 
something like the figure below. (I believe I used nr=4. You are encouraged to 
experiment with propensities and run numbers.) 



3G 




time 

Figure 8.4. Dimer Self-control. 
Your work will be graded as follows: 

6 pts for header CONTAINING detailed USAGE 
6 pts for further comments in code 
4 pts for indentation 
10 pts for correct gilldriver 

10 pts for correct computation of reaction number 
8 pts for correct coding of update 



6 pts for two plots corresponding to the free dimer 

production, and to its self-regulated production. 
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9. Deterministic Chemical Kinetics 



The stochastic view is a bare-handed means of following the dynamics of indi- 
vidual molecules. For large reacting systems this approach may be prohibitively 
expensive and/or unnecessarily precise. In such cases one may rely on more 
coarse grained tools, e.g., the 

Law of Mass Action: The rate of a chemical reaction is directly proportional 
to the product of the effective concentrations of each participating reactant. 

We will come to understand this law by its application. 

Example 1. Suppose R\ is the bi-reactant system 

X\ + X2 — > X3 

that proceeds with propensity c\ in a region of volume V. It follows that there 
are X1X2 distinct combinations of reactants inside V, and so X\X<iC\dt gives 
the probability that R\ will occur somewhere inside V in the next dt units of 
time. If we now average over nr runs, 

mean (X\X2C\) = c\ mean (X1X2) 

is the average rate at which R\ reactions are occurring inside V. The average 
reaction rate per unit volume is therefore 

mean {XiX^Ci/V = Fcimean (X1X2) 

in terms of the molecular concentrations 

X 3 EE Xj/V. 

Now, the reaction rate constant, ki, is defined to be this average reaction rate 
per unit volume divided by the product of the average concentrations of the 
reactants 

Fcimean {x\X2) 
1 mean (xi)mean (#2) 
Finally, on equating the average of a product with the product of the averages, 
i.e., 

mean {x\X2) = mean (xi)mean (#2) 

we arrive at 

fci = Vcl (9.1) 



: J »8 



In this case the law of mass action dictates that 

x[(t) = x' 2 (t) = — k\X\{t)x2{t) and x' 3 (t) = k\Xi{t)x2{t) . 
Example 2. Suppose R2 is the mono-reactant system 

that proceeds with propensity C2 in a region of volume V. In this case there 
are X\(X\ — l)/2 possible pairs and if X\ is large then this number is relatively 
close to Xf/2. As such 

k 2 = Vc 2 /2 

and mass action dictates that 

x'i{t) = —2k,2x\{t) and x' 3 (t) = k,2x\{t). 

Example 3. We return to (8.3) 

-^00 + y —> %Y (propensity c\) 
2Y — > Z (propensity c 2 ) 

and, with k\ = c\V and k2 = C2V/2 find 

y'(t) = 2k\x O0 y{t) - (ki (t) + 2k 2 y 2 (t)) 
gains — losses 



(9.2) 



If we denote the initial concentration y(0) = yO we may solve the ordinary 
differential equation (9.2) "by hand" (via calculus or 

dsolve('Dy = klx*y-2 *k2 *y " 2 ' , ' y ( ) =y0 ' ) 

in MATLAB) and find 

/,s yohx^ 

y\t) = 



2y k 2 - (2y k 2 - hx^) exp(-/cix 00 i) 

This function is easy to graph and so compare with its stochastic cousin. As- 
suming V = 1 we arrive at the figure below on running gilllo.m 
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time 



Figure 9.1. Stochastic vs. Deterministic approach to (8.3). 

We see that the deterministic model does an excellent job of tracking the 
stochastic trajectory. 
Example 4. The Lotka system 



Y x + Y 2 
Y 2 



2Y\ propensity c\ 
2Y 2 propensity c 2 
Z propensityc3 



CjV : 



yields, with kj 

y[ = hxooVi - k 2 yiy 2 

y' 2 = hym - hy2 

The solution of this system is not presentable in "closed form" and so we turn to 
approximate, or numerical, means. Matlab sports an entire suite of programs 
for solving differential equations. The default program is ode2 3. We invoke it 
in our Lotka 2-way code and arrive at the figure below. As above blue is the 
Gillespie trajectory and red is the ode trajectory. 

Lotka Two Ways 




500 1000 1500 2000 2500 

number of y 1 molecules 



Figure 9.2. Stochastic vs. Deterministic approach to (8.4). 
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Example 5. The Brusselator system 



^00,1 

2Y 1 + Y 2 
Y 1 



Y\ propensity c\ 
Y2 + Z\ propensity C2 
3Y\ propensity C3 
Z2 propensity C4 



becomes 



(9.3) 



y[ = ci^oo,! - 02X00,22/1 + (3 - 2)(c 3 /2)2/i2/ 2 - cm 
y' 2 = €2X00,22/1 - {c 3 /2)yly 2 

This associated ode code produces the figures below. You may wish to zoom in 
to better see the transitions that each species undergoes, or better yet, take a 
look at the (2/1,2/2) plane. 
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Figure 9.3. Stochastic vs. Deterministic approach to (9.3). 

Deterministic vs. Stochastic Operon Simulation 

We continue last week's investigation of a self-governing gene. Recall the 8 
reactions in (8.5) and (8.6) involved the 6 species 

X =[RPr RPr P D I] 

and propensities c(l) throught c(8). 

Your first task is to write down the associated 6 ordinary differential equations 
(in the order they are listed in x above). I will give you the first one, 



Jb 1 



■C\X\X2 + C 2 X 3 + C3X3 



Your next task is to augment last week's code so that it solves (via ode2 3) 
this ode system and compares the mean (stochastic) dimer level to the ode 
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prediction. For example, using the initial counts and propensities 

xO = [10 1 0] and c = [2 1 4 2 0.5 0.05 .01 .005] 
and nr = 10, we arrive at the figure below 



250 
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20 40 60 80 100 
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Figure 9.4. Average Stochastic vs. Deterministic approach to the self- 
regulating gene. 

This result is promising, in that the dimer predictions appear to coincide. It 
would be nice to dig a little deeper and ask how the other reactants are fairing. 
At the stochastic level we know that each reactant may only take integer jumps 
and that R, Pr, RPr and I in fact jump back and forth between and 1. As 
taking ensemble means introduces intermediate values we instead compare our 
ode solution to single runs of Gillespie. The figure below was obtained with the 
same initial values and propensities as above. 
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Figure 9.5. Solo Stochastic vs. Deterministic approach to the self-regulating 
gene. 
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This figures cries out for explication. It shows early switching in R, Pr and 
RPr, associated with stable early growth in D, despite the high volatility in P. 
And finally, once D achieves some critical level we see I switch on and stay on. 

You will accomplish these two figures by changing last week's gilldriver to 
a function called mca2driver that sets sets tfin, tine, nr, xO, and c (first with 
nr=4 or so), calls mygill nr times and plots the mean and standard deviation 
like above, then calls ode23 which in turn calls 

function dx = mcaode (t,x, c) 

and then plots the new dimer count as in our first figure above. Your driver 
then sets nr = 1 and constructs, via subplot (2, 3, j) , the lovely six panel 
figure above. In order to "see" through the noise I have limited the width of 
this window relative to the tfin used above. You too can do this with xlim ( [0 
30] ). 

Your work will be graded as follows: 

6 pts for header CONTAINING detailed USAGE 

header of mcaode should give full system in detail 
6 pts for further comments in code 
4 pts for indentation 
14 pts for correct mca2driver 
10 pts for correct mcaode 

10 pts for two plots similar to those above 
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10. Modeling of Fiber Networks 



We consider the case illustrated in Figure 10.1. The bold solid line is a fiber 
in its reference state. When we subject its two ends to the two forces, (/i,./2) 
and (/3,/4) the respective ends are displaced by (xi,X2) and (#3, £4). 




Figure 10.1. The reference (solid) and deformed (dashed) fiber. 

Our goal is to build a theory that predicts x from knowledge of /. The 
first step is to quantify the associated elongation, e = I — L, where L is the 
undeformed length and £ is the deformed length. With respect to Figure 10.1, 
we suppose that the lower left node of the undeformed fiber sits at (0, 0) in the 
Cartesian plane, while its upper right node resides at (Lcos#, Lsin#). Following 
Euclid, we write 

I = \J (L cos 9 + £3 — xi) 2 + (L sin 6 + X4 — X2) 2 
= \/ L 2 + 2L{(x3 — x\) cos# + (X4 — X2) sin9} + (x% — xi) 2 + {x± — X2) 2 
= L\Jl + 2{(x3 — x\) cos 6 + (1C4 — X2) sin 9}/L + {(x^ — x\) 2 + (x~4 — X2) 2 }/ ' L 2 . 



To help us here we invoke MacLaurin, y/1 +t = 1 + 1/2 + 0(t 2 ) for small t, and 
write 

I = L + (xs — x\) cosO + (x4 — X2) sin# + 0((xj — Xj + 2) 2 / L). 

Hence, assuming that (xj — Xj + 2) 2 is small (for both j ' = 1 and 2) compared with 
L, we find 

^3 — Xi) cos ^ + (x4 — X2) sin ^. (1) 
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In the future we shall write = instead of It will be understood that we are 
working under the hypothesis that the end displacements are small in comparison 
to the undeformed length. This approximation can be interpreted geometrically, 
or visually, to mean simply that each fiber tends to deform along its original 
direction. 

We assume that our fiber is Hookean, in the sense that its restoring force, y, 
is proportional to its elongation. More precisely, we presume that 

Ea 



where E denotes the fiber's Young's modulus, a denotes its cross sectional area, 
and L denotes its reference length. This y, positive when the bar is stretched 
and negative when compressed, acts along the reference direction, 9, in balance 
with the applied load /. More precisely, at the lower node 

y cos 9 + fi = and y sin 9 + = (3) 

and at the upper node 

y cos(7r + 9) + /a = and y sin(7r + 9) + f± = 

or 

-y cos(#) + / 3 = and - y sin{9) + / 4 = (4) 

Finally, we need only substitute our expression for y in terms of e and e in terms 
of x and recognize that the 4 equations in (3) and (4) (hopefully) determine x 
from /. More precisely, with k = Ea/L, we find 

k{(xz — x\) cos 9 + (#4 — X2) sin9} cos 9 = —f\ 
^{(^3 — x\) cos 9 + (x4 — X2) sin9} sin# = —ji 
k{{xs — x\) cos 9 + [xa — X2) sin 9} cos 9 = f% 
k{ (^3 — x\) cos 9 + (x4 — X2) sin 9} sin 9 = f± 

Let us consider a few concrete examples. If 

k = 1, 9 = 0, h = h = and h = -/ 3 

the above system of four equations deflates to 

xz - xi = h 

which indeed determines x\ and X3 up to an arbitrary rigid motion, stemming 
from the fact that our fiber is a floater. 
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Example 1. We nail things to a foundation, and add a fiber and arrive at the 
tent below. 




Figure 10.2. A simple tent. 
We first compute the two elongations 

ei = x\ cos(#i) + x 2 sin(#i) 
e 2 = x\ cos(# 2 ) + x 2 sin(# 2 ) 

we suppose the fibers have stiffnesses k\ and k 2 and so 

yi = k x ei 
2/2 = k 2 e 2 

while force balance at the only free node yields 

-yi cos(0i) - y 2 cos{6 2 ) + /i = 
-yi sin(^i) - y 2 sin(# 2 ) + h = 

Assuming B\ = 7r/4 and 9 2 = 37r/4, we find 

ei = (xi + x 2 )/v2 
e 2 = (x 2 - xi)/v2 

and 

0i - 2/2 )/v^ = /i 
(2/1 + 2/2)/V2 = / 2 

and so x must obey 

(fci + /c 2 )xi/2 + (A?! - k 2 )x 2 /2 = /1 
(fci - fc 2 )zi/2 + 0i + ^2)^1/2 = / 2 
In the case of bars of equal stiffness we find the simple answer that 

xi = fi/h and x 2 = f 2 /k 2 . 
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For unequal stiffnesses we must solve our 2 linear equations simultaneously. 
Matlab knows Gaussian Elimination. 

We are interested in understanding big nets (say m fibers meeting at n joints) 
and so we step back and realize that our model was constructed in three easy 
pieces. 

The fiber elongations are linear combinations of their end displacements, 

6 == Ax , 

where A is m-by-2n, is called the node-edge adjacency matrix, and encodes the 
'geometry' of the fiber net. 

Each fiber restoring force is proportional to its elongation, 

V = K e, 

where K is m-by-m and diagonal and encodes the 'physics' of the fiber net. 
The restoring forces balance the applied forces at each node, 

A T y = f 

where A T is the transpose (exchange rows for columns) of A. 
When these steps are combined, we arrive at the linear system 

A T KAx = f. 

For the net of Figure 10.2, we have 



A 



'cos 6*i sin 9\ 
cos 6>2 sin 6*2 



K 



'Jfei o N 
k 2 , 



A 



T 



'cos 9\ cos 62 
sin 61 sin 62 



We have coded this tent example here. 



Example 2. We now build the adjacency matrix for the trailer below. 

I 3 ' 2 6 




Figure 10.3. A trailer. 
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We have numbered the 7 fibers and the 3 nodes. We shall adopt the conven- 
tion that the horizontal and vertical displacements of node j are x<ij-\ and X2j 
respectively. With the fiber angles, 

1 = 4 = 7 = tt/2, 6 2 = 6 5 = 7t/4, and 6 3 = 6 6 = 0, 

the associated elongations are 

ei = x 2 

e2 = {%3 + X4)/V2 

e3 = X'3 ~ Xl 
64 = £4 

es = (x 5 + x 6 )/V2 
e 6 = x 5 - Xi 
e 7 = x 6 . 

which we translate, 1 row, i.e., one fiber, at a time. 
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We have coded this trailer example here. We next animate this scene by applying 
an incremental force in our Moving Trailer Code 



Project: Plaza Castilla Elastica 

You will load, solve and animate a class of cantilevers like that exhibited 
below. This one is 4 stories tall and so has 4x2=8 nodes (and so 2x8=16 degrees 
of freedom) and 4x4=16 fibers. The odd numbered fibers make an angle of 7r/4 
with the horizontal. You will write a function 

function cantilever (Ea, F,nof ) 

where Ea is the 4xnos-by-l vector of modulus area products, F is the strength 
of the downward force at the upper right tip, and nof is the number of frames 
in the resulting collage. We have denoted the number of stories by nos. 
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/ 10 8 




5 10 15 

degree of freedom 



Figure 10.4. A cantilever and its adjacency matrix. 

If you follow my ordering of the nodes and fibers then your adjacency matrix 
will exhibit considerable structure. This structure is revealed through Matlab's 
spy command as in the figure at the upper right. 

Your code must accomodate an arbitrary number of stories and it must pro- 
duce 2 figures. The first figure should come from spy (like that at the upper 
right) while the other should be a likeness of the Duchampian collage below. 
This figure is obtained by scaling the vertical load by frame number, as in our 
moving trailer example. I used one call to line to draw the net at each load 
and I used fill to draw the base. 




Figure 10.5. A descending cantilever. 
Your work will be graded as follows: 
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6 pts for header CONTAINING detailed USAGE 

8 pts for further comments in code 

4 pts for indentation 

8 pts for correct undeformed coordinates 

8 pts for correct adjacency matrix 

8 pts for correct plotting of deformed cantilever 

4 pts for labeled plot of adjacency matrix for nos = 6 

4 pts for plot of Duchampian cantilever for nos = 6 
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11. Gaussian Elimination 



Looking back at our fiber code we see that all of the 'math' is hidden under 
the innocuous backslash - namely x=s\f . We strive in this subsection to reveal 
the hidden work by examining systems of equations of increasing complexity. 

The easiest systems to solve are the uncoupled, or diagonal systems, i.e., 
systems of the form Sx = f where only the diagonal elements of S are nonzero. 
If S is n-by-n we simply write 

for j=l : n 

x(j) = f (j)/S(j, j); 

end 

or, even simpler, 
x = f . /diag ( S ) 

Of course this procedure breaks down if S has a zero on its diagonal. A diagonal 
matrix with one or more zeros on its diagonal is said to be singular . In such 
a case, Sx=f is not solvable unless f has a corresponding zero. Although, even 
in this case we note that x is left undetermined. For example 




determines x(l) but not x(2). 

After diagonal matrices the next most easily handled type is the class of 
triangular matrices. These are matrices all of whose nonzeros lie either on and 
above (upper) or below (lower). In the following upper triangular system 




we start at the bottom and notice that x(2) = 6/3 = 2. The first equation then 
reads 4r(l) + 2a?(2) = 4 or 4x(l) = 4 - 2x(2) or x(l) = (4 - 2a? (2)) /4 = 0. This 
procedure is commonly known as backsubstition. With respect to coding, it 
requires only slightly more effort than in the diagonal case. 

x = zeros (n, 1) ; 
x(n) = f (n) /S (n,n) ; 
for j=n-l : -1 : 1 
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tmp = ; 
for k= j+1 : n 

tmp = tmp + S(j,k)*x(k); 

end 

x(j) = (f (j) - tmp)/S(j, j); 

end 

As in the diagonal case we find ourselves dividing by each diagonal element of S . 
Hence, a triangular matrix with one or more zeros on its diagonal is said to be 
singular . Gaussian elimination brutally seeks to reduce a given linear system 
to a triangular linear system via repeated application of two elementary row 
operations. The first is Row Swapping, for example, 

S([j k], :) = S([k j],:) 

and the second is Row Mixing, for example 

S(j,:) =S(j,:) + magicnumber *S (k, : ) 

In each column we swap in order to bubble up a big pivot then we repeatedly 
mix until the pivot has eliminated every element below the diagonal. The 
systematic use of these two utilities is diagrammed in the following pseudo code: 

function x = gauss (S,f) 
n = length ( f ) ; 

S = [S | f] Augment S with f 

for k=l:n-l k counts columns 

r = row number, larger than or equal to k, 

with largest value (in magnitude) in column k 

if this largest value is really small then warn the user 

swap row r and row k 

for j=k+l : n 

mix row k into row j in order to eliminate S(j,k) 

end 

end 

if S(n,n) is really small then warn the user 

strip off the changed f, i.e., copy column n+1 of S onto f 

x = trisolve (S, f ) 

return 

Here is a small example where S and f are: 
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2 



4 



2 



2 
3 



1 



2 



4 



and 



4 



2 



1 



1 



Here is S augmented with f: 

2 4 2 2 
12 4 3 
4 2 11 

rows 1 and 3 are swapped: 

4 2 11 
12 4 3 
2 4 2 2 

our first pivot is 4 and elimination occurs in column 1: 

4 2 1 1 

3/2 15/4 11/4 

3 3/2 3/2 

rows 2 and 3 are swapped: 

4 2 1 1 

3 3/2 3/2 

3/2 15/4 11/4 

our second pivot is 3 and elimination occurs in column 2: 

4 2 1 1 

3 3/2 3/2 

3 2 

our third and final pivot is also 3, our matrix is now triangular, and so we may 
trisolve: 



4*x(l) + 2*(l/6) + l*(2/3) =1 so x(l) = 

In the case that you encounter a small (smaller than eps) pivot your warning 
message might look something like that found in my singular diary. 



x(3) = 2/3 

3*x (2) + (3/2) * (2/3) = 3/2 



so x(2) = 1/6 



53 



Project: Fiber Nets and Gaussian Elimination 

Nonzero adjacencies, cx = 6 cy : 
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Figure 11.1. A cantilever and its adjacency matrix. 

You will load, solve and animate a class of cantilevers like that exhibited 
above. This one is 6 cells wide and 4 cells tall. It has 102 fibers and 60 degrees 
of freedom. If the fibers and nodes are numbered as indicated then one arrives 
at an adjacency matrix like that spied above. 

Regardless of the number of cells, this cantilever has length 2 and height 1. 
As a result, the angles of the cross bars will depend upon your choice of cx (the 
number of horizontal cells) and cy (the number of vertical cells). NOTE: The 
diagonal bars cross over one another without calling for a node. 

To make sure we are all on the right track let us build the adjacency matrix 
by hand when cx = cy = 2. In this case, see fig. we have 18 fibers connected at 
6 nodes and so the net enjoys 12 degrees of freedom. 




Figure 11.2. The cantilever with cx = cy = 2. 
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We assess their elongations via repeated reference to our general formula, 
equation (1) on page 44. 

To begin, the elongation of the six horizontal fibers obey 

ei = a?i, e 5 = x 3 - xi, e 9 = a? 5 , 
ei3 = — X5, ei7 = xg, eis = % — £9, 

while the elongations of the vertical fibers obey 

e4 = xq — X2, eg = xg — X4, ei2 = a?io — xq, e\§ = xyi — xg. 

Regarding the diagonal fibers we must first compute the associated angles. The 
length of fiber one is (2/cx) = 1 while the length of fiber 4 is 1/cy = 1/2 and so 

cos(#i) = 2/v / 5 and sin(^) = 

With s — l/\fb equation (1) now provides 

e 2 = 2sx 5 + sxq, e 6 = 2s(x 7 - x\) + s(xg - x 2 ), 
eio = 2sxq + sxio, ei4 = 2s(xn - x 5 ) + s(xu - x 6 ). 

Finally, as 62 = tt — 9i we find that 

cos(# 2 ) = -2/Vb and sin(# 2 ) = 1/Vb 

and so, 

e3 = — 2s(0 - xi) + s(0 - x 2 ), e 7 = -2s(x 5 - x 3 ) + s(x 6 - £4), 
en = -2s(0 - x 5 ) + s(0 - x 6 ), ei 5 = -2s(a:9 - x 7 ) + s(rci - x 8 ). 
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On gathering all 18 elongations we arrive at the adjacency matrix 



A = 



1 
































0\ 














2s 


s 




















2s 


— s 



































-1 











1 




















-1 





1 





























-2s 


— s 














2s 


s 




















2s 


— s 


-2s 


s 





























-1 











1 


























1 















































2s 


s 




















2s 


— s 



































-1 











1 




















-1 





1 





























-2s 


— s 














2s 


s 




















2s 


— s 


-2s 


s 





























-1 











1 


























1 



































-1 





1 


0/ 



I recommend that you test your code against these values. 
You will adapt and extend last week's function to 

function cantilever 2 (cx, cy, Ea, F, not ) 

where cx and cy dictate the geometry of the net, Ea is the vector of modulus 
area products, F is the strength of the downward force at the center of the 
right face (node 18 in the figure above), and nof is the number of frames in the 
resulting movie. NOTE: The right face has no center unless cy is even. Your 
code may presume that the user will only call it with an even cy. 

In addition to encoding a new net this week's code must also replace last 
week's miraculous x=s\f with 

x = gauss ( S , f ) 

where gauss is a subfunction that you must write and append to cant i lever 2 
Recall that this subfunction has been sketched for you on the lecture page. In 
particular, it itself calls the subfunction trisolve. You should of course test 
your gauss against backslash on a number of random matrices before tucking it 
into cantilever2. 



56 



Figure 11.3. The bending of the cantilever. 

Your code must accomodate an arbitrary number of cells and it must produce 
a figure via spy (A) like that above, and a Duchampian trace like that at right, 
of a net undergoing an incremental force at the center of its right face. (F = 
0.05, nof = 5). 

Your work will be graded as follows: 

6 pts for header CONTAINING detailed USAGE 
6 pts for further comments in code 
4 pts for indentation 

4 pts for correct undeformed coordinates 
6 pts for correct adjacency matrix 

6 pts for correct plotting of deformed cantilever 
4 pts for correct row swapping in gauss 
6 pts for correct row mixing in gauss 

4 pts for labeled plot of adjacency matrix, cx=4, cy=6 
4 pts for 5 frame plot of moving cantilever, cx=4, cy=6 
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12. Optimal Design of Fiber Nets 

Biological fiber nets have (likely) evolved to best fit their niche. How one 
measures 'fitness' and how one identifies nature's 'design variables' are each 
matters of considerable debate. We will focus on the simplest such measure, 
namely compliance, or work done by the load, and suppose that nature may 
tune the radii of her fibers. Recalling that work is merely the product of force 
and distance we express compliance as 

C = x T f (12.1) 

where x is the solution to 

A T K{a)Ax = f (12.2) 

where we have stressed the dependence of K, the elemental stiffness matrix, upon 
a, the vector of fiber cross sectional areas. We will presume throughout that 
each fiber has the same Young's modulus. Our interest then is in minimizing the 
C of (12.1) subject to x obeying the equilibrium equation (12.2). A moment's 
reflection will confirm that the compliance can be reduced to zero by simply 
choosing infinitely fat fibers. We preclude this possibility by constraining the 
volume of the net, via 

L T a = V (12.3) 

where L is the column vector of fiber lengths. 

Let us solve a few small design problems by hand before opening the big 
toolbox. In particular, the displacement of the loaded symmetric (9 = 7r/4) tent 
is 

fx(l)\ = s fa(l) + a(2) a(2) - a(l)\ ff(l)\ 
\x(2)J fl (l)a(2) \a{2) - a(l) a(l) + a(2) J \f(2)J 

where s 2 = 1/2 as usual. And so 

c = 5 (/(l)-/(2)) 2 + S (/(l) + /(2)) 2 



a(2) a(l) 
If we now invoke the volume constraint 

a(l) + a(2) = sV 

we find that 

5 (/(l)-/(2)) 2 5 (/(l) + /(2)) 2 
sV-a(l) a(l) 



58 



and so we have reduced ourselves to a one-dimensional minimization problem. 
Let us consider a pair of special loadings. 

Design against falling debris: In this case / = [0; —1] and so 

s s _ s 2 V 

~ a(T) + sV - a(l) ~ a(l)(sV - a(l)) 

which attains its (positive) minimum at 

a(l) = a (2) = sV/2. 

This tells us that the two tent legs should be equally fat to best sustain a vertical 
load. 

Design against bears: In this case / = [1; 1] and so 

4s 

C = 



a(l) 

which approaches its minimum as a(l) approaches infinity. This of course leads 
to the absurdity that a(2) approaches negative infinity. If we impose the ad- 
ditional constraints that a(l) and a (2) each be nonnegative then we arrive the 
optimal design 

a(l) = sV, a(2) = 0. 

This corresponds to the fact that, as the load is perpendicular to the tent's right 
leg, the right leg does no work. In reality, in order to account for variations in 
the loading, one typically sets a nonzero lower bound on each a(J), and for cost 
or construction contingencies one likewise also imposes a finite upper bound on 
each a(j). In other words, we ask that 

alo < a(j) < ahi for each j. (12.4) 




Figure 12.1. Best design against bears. 
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We now arrive at the General Case 



minx / 

a 

subject to 

A T K(a)Ax = /, L T a = V and alo < a(j) < ahi for each j. 

Matlab has a large toolbox for solving such optimization problems, type 
help optim. The routine most pertinent to our problem is fmincon. Its 
application to our simple tent is coded here. Running this code with / = [.5; — 1] 
reveals the figure at right. 

As we are now moving into larger nets we should attempt to give fmincon 
a hand, fmincon searches for a minimum of the constrained compliance by 
searching for a critical point, that is an a for which the gradient, i.e., the vector 
of derivatives of the compliance with respect to each a(j), is zero. Hence, to 
accelerate fmincon we need only pass along information on the gradient of the 
compliance. The gradient we need is 



VC(a) 



f dC{a) dC{a) dC{a) 
\ da\ dd2 da m 



T 



Let's compute these one at a time, e.g., 



dC(a) j, T dx(a) 



dai da\ 

To compute the latter we differentiate the equilibrium equation 

A T K{a)Ax{a) = f 

with respect to a\ and find, via the product rule and the fact that A and / are 
independent of a, 



We solve this for 



A 1 —^Ax(a) + A 1 K(a)A—^- = 0. 

OCL\ OCLi 
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and so arrive at the desired 

„ T dx(a) T 



r 



da\ 



= -x T (a)(A T K(a)A)(A T K(a)A)- 1 A rdK{a] 



Ax(a) 

(Xjhx 




fi\ e 2 



V 



o o 







da i 



/ ei \ 

&2 



Ax(a) 



= -e?(a)/L!. 

Arguing in a similar fashion for the remaining partial derivatives we find 



VC(a) 



[e\{a)/Li e 2 2 {a)/L 2 



e m (a)/L m 



Do you see how minimizing the compliance is synonymous with minimizing the 
square of each elongation? 

Project: Optimal Design of Fiber Nets 

We turn to fmincon to design the strongest possible cantilever of volume V 
with the option of freezing the size of the horizontal fibers. 

Compliance = 0.86884 




Figure 12.2. The best cantilever in its class. 

You will write a function 

function cantilever 3 (cx, cy, F, V, nohor ) 
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that builds last week's adjacency matrix, A, vector of fiber lengths, L, and force 
vector, f . This f should be all zeros except for the value F in the component 
corresponding to vertical force at the center of the right face. (There is no 
center if cy is not even). Your cantilever3 will also generate an initial guess for 
the fiber cross-sectional areas, aO, that jibes with the volume constraint, 1/ *a0 
= V. Your bestcant will then set lower and upper bounds on the permissible a 
and then call fmincon to find the best a. fmincon will in turn call your 
compcant function, where the compliance of the fiber net is computed. Each 
of these steps is a straighforward generalization of our best tent code. 

You will inform fmincon that you are supplying the gradient of your objec- 
tive function by switching on one of its options and by coding compcant to 
return both the compliance (1-by-l) and its gradient (m-by-1). In particular, 
you must follow this protocol 

function [comp,grad] = compcant (a, A, L, f) 

as in the simple example. 

Once fmincon returns the best a it should be plotted, as above. For visual 
purposes you may wish to use the fat-factor-5 of btent. 



Compliance = 1.0433 




Figure 12.3. The best cantilever in the class of cantilevers without horizontal 
stiff eners. 

Finally, you will write a simple driver 
function cant3drive 
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that sets cx, cy, F and V as well as specifies whether or not all of the horizontal 
fiber areas should be held at their LB values (if nohor=l then UB=LB=0.001 
for each horizontal fiber). With nohor=0 you should achieve the figure above. 
With nohor=l you should achieve the figure at right. 

Your work will be graded as follows: 

6 pts for header CONTAINING detailed USAGE 

8 pts for further comments in code 

4 pts for indentation 

8 pts for correct cantilever3 

8 pts for correct compcant with gradient 

with cx = 4, cy = 6, V = 20 and F = -0.5 

8 pts for labeled plot of best net, nohor=0 
8 pts for labeled plot of best net, nohor=l 
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13. Multi-Layer Perceptrons 

As our next example we shall study networks of simple minded neurons. At- 
tractors on these nets will look and behave much like the attractors for our gene 
nets. Our interest however is not the behavior of a single net but rather a net's 
capacity to learn in both supervised and unsupervised settings. 

This is not so far from your work on optimal fiber nets. In that case you 
began with a random set of fiber diameters and a given load and the net slowly 
learned how to redistribute its volume in such a way as to minimize the work 
done by the load. The role of teacher was played by fmincon. 

In this setting we shall teach a neural network to respond in a desired way 
when faced with certain stimuli. 

To learn is to modify one's synaptic weights. 

Let us start with a simple example. One neuron with two inputs. The neuron 
scales and sums its inputs and fires if the sum exceeds a given threshold. If the 
inputs and output are binary and the threshold for firing is 0.5 then 

out = round (w ( 1 ) *in ( 1 ) + w(2)*in(2)) 

is our neural net. The challenge is to choose wl and w2 in accordance with 
a particular task. For example, can we teach this neuron to perform the and 
operation? More precisely, we search for wl and w2 such that 

round(w(l)*0 + w(2)*0) = 

round(w(l)*0 + w(2)*l) = 

round(w(l)*l + w(2)*0) = 

round(w(l)*l + w(2)*l) = 1 

Here is one possible approach. You see that w learns to associate each stimulus 
with its preferred target by following the negative gradient of the error. 

This example is of course a few hundred billion cells short of a full brain and 
is perhaps a little abrupt with respect to threshold. Let us consider bigger more 
gentle networks. We will consider networks with 3 layers comprised of 

n inputs 

h hidden cells, and 
m outputs 

The h-by-n weighting matrix from inputs to hidden cells will be called V and 
the m-by-h weighting matrix from hidden cells to outputs will be called W. We 
offer a concrete example in the figure below. 
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Figure 13.1. A three-layer Multiperceptron. 
The hard threshold at each cell will be replaced by the soft sigmoid 

s (x) = 1 . / (1 + exp (0 . 5-x) ) 
With this notation, if p is an input pattern then the output will be 
o = s (Ws (Vp) ) 

This representation is dangerously brief - please make sure you understand it 
before proceeding. 

We will typically have N input patterns that we lay into the columns of the 
n-by-N matrix 

p ( : , 1) , p ( : , 2) , p ( : ,N) 

that we would like to pair with N targets that comprise the columns of the 
m-by-N matrix 

t (: , 1) , t (: ,2) , . . . , t (: ,N) . 

We do this by choosing V and W to minimize the misfit 

1 N 

i=l 

1 N 

= 2 E^W^'O)) - i)) T (s(Ws(Vp(:M - t(:, i)) 

i=l 
1 N m 

= 2 E EW^'- : WVp(:, »))) - i)) 2 . 

i=l j=l 

We minimize this sum of squares by guiding V and W along the gradient of E. 
That is, we update V and W according to 

W = W — rgrdud w E 

V = V — rgrad^E" 
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for some fixed learning rate, r. These derivatives are easiest to see if there is 
but one input pattern (N=l) and a single output cell (m=l), for then the misfit 
is simply 

E{V, W) = (l/2)(s(Ws(Vp)) - tf 
and so the chain rule reveals the gradient, or vector of partial derivatives 

grad^£(V, W) = [dE/dWi dE/dW 2 ■ ■ ■ dE/dW h ] 

takes the form 

grad_W E(V,W) = ( s (Ws (Vp) ) -t ) * s' (Ws (Vp) ) * s (Vp) ~T 

lxl lxl lxh 

And similarly, the matrix of partial derivatives, 

fdE/dV^i dE/dV^ ■■■ dE/dV 1>n \ 
dE/dV 2 ,i dE/dV 2 , 2 ■■■ dE/dV 2 , r . 



grad v E(V, W) 



\dE/dV hA dE/dV hj2 ■■■ dE/dV h J 



takes the form 



grad_V E (V, W) = (s (Ws (Vp) ) -t) *s' (Ws (Vp) ) * (W~T. *s' (Vp) ) *p~T 

lxl lxl hxl lxn 

These expressions can be simplified a bit upon observing that our sigmoid obeys 

s' (x) = s (x) (1-s (x) ) 

and so, setting 

q = s (Vp) and o = s (Wq) brings 

grad_W E = (o-t) o (l-o) q"T and 

grad_V E = (o-t)o(l-o) (W^T . *q. * ( 1-q) ) p^T 

We have coded this training here. We test the net it generates using nnxor, in 
this diary 
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Figure 13.2. A multiperceptron that performs exclusive OR. 



>> [V,W] = xortrain([5 5;1 1],[10 -15 ] , 10000 , . 25 ) 
V = 6.0321 6.0311 
1.0216 1.0066 
W = 13 .4652 -17 . 7826 

>> nnxor ( [0 0] ' , V,W) 
ans = 0.1062 

>> nnxor ( [0 1] ' , V,W) 
ans = 0.8600 

>> nnxor ( [1 0] ' , V,W) 
ans = 0.8524 

>> nnxor ( [1 1] ' , V,W) 
ans = 0.1615 

Do you see that our pupil has stumbled upon the Boolean Identity 

XOR(a,b) = OR(a,b) - AND ( a , b ) 

You see from the code that we have proceeded as if there was only one input 
pattern by simply choosing one at random each time. The same ploy may be 
used in the case of multiple outputs, as, e.g., in this week's assignment. 

Project: - DNA OCR via Neural Nets 

Train a 3 layer neural net to recognize digitized versions of the letters A, C, 
G and T. Your net should have 25 inputs, 2 outputs and 25 "hidden" cells. The 
inputs will be binary and correspond to the LED-like figure below. 
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2 


3 


4 


5 


6 


7 


a 


9 


10 


11 


12 


13 




15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 



Figure 13.3. The input pattern [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]. 
For your targets, please use 

T = [0 0], A = [0 1 ] , C=[10], and G = [1 1]. 

Write a function 

[V,W] = ocrtrain (V0 , W0 , maxiter , rate) 

that returns the two weight matrices after successful training at the given rate. 
Drive this function with a function called ocrdrive that 

sets maxiter=5000 and rate=0.1, 

sets initial weights V0 and W0 via randn, 

sends all this info to ocrtrain and receives 

the trained V and W, 
displays V and W as below via imagesc, colorbar, 

hist, reshape, and subplot 
applies these learned V and W to a 12-by-25 

bitstream whose 1st four rows are 
[1 11111000010000100001111 1; 
111111000010011100011111 1; 
111111000111111100011000 1; 
111110010000100001000010 0] 
to get the next 4 four rows flip 1 bit in each of the above, 
to get the final 4 four rows flip 2 bits in each of the above. 
Displays these 12 letters as below via, subplot, 

imagesc, reshape and colormap ( 1-gray ) , 
Finally, display in the command window the sequence of letters 
that V and W determined were represented in the bitstream, 
e.g., seq = CGATCGATCGAA 

Hint The grad formulas established above work nice for scalar outputs. But 
here your outputs have two components, in particular your W has two rows. 
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What to do? Well, on each iteration just work on updating a single row of W. 
How should you choose which row? Flip a rand. 2nd Hint: To get exactly my 
V and W we should start at the same random VO and WO. To do that please 
reinitialize randn before each use via randn ( ' state' , ) . 

Trained V Distribution of trained V values 

>:■■.■_.::.!,■ 

Trained W 




Figure 13.4. Representations of the trained weight matrices. 
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CERT 



Figure 13.5. Progressively noisier input patterns. 
Your work will be graded as follows 

8 pts for headers CONTAINING detailed USAGE 

8 pts for further comments in code 

4 pts for indentation 

10 pts for correct ocrtrain 

8 pts for correct ocrdrive 

4 pts for V and W titled imagesc plots 

4 pts for V and W titled hist plots 

4 pts for display of final letters to command window 
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14. Hopfield Nets 

Hopfield nets are comprised of N hard-threshold nodes with all-to-all sym- 
metric coupling where on = 1 and off = -1. Let us start with the net below 




Figure 14.1. A four node Hopfield net. 

There are 4 nodes and so we must build a 4-by-4 weight matrix W. The bidi- 
rectional arrows imply equal weights in each direction, i.e., 

W(iJ) = W(j,i) 

If s is the current state of the net then the next state is 

ns = sign2(VFs) 

where 



sign2(rc) = 



1 if x > 
-1 if x < 



is our hard threshold function. This net can be trained to remember an input 
pattern p by setting the weights to 



T 



W = pp 

then p and its mirror image —p will be attractors. That is, they will satisfy 
s = sign2(M /r s). To see this note that (as p T p > 0) 



Wp = pp T p = pip 1 p) = (p 1 p)p and so sign2(Wp) = sign2(jp) = p. 



T 
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In fact, for each initial state, s, the very next state will be either =bp or 
-ones (N, 1 ) (when p T s = 0), for 

Ws = pp T s = p(p T s) = (p T s)p and so sign2(iys) = sign2(j> T s)sign2(j>), 

unless p is orthogonal to s, i.e., p T s = 0, in which case 

sign2(VFs) = — ones(N, 1). 

If this full off vector is itself orthogonal to p then it and ±]9 will be the only 
attractors. To get a feel for this I encourage you to dissect and exercise hop on 
several 4-by-l choices. 

All of this generalizes nicely to multiple training patterns. In fact, if p(:, 1) 
and £>(:, 2) are two such patterns we set 

P= [p(:,l) p(:,2)] and W = PP T 



Arguing as above, we find 



Ws = PP 1 s = l))p(:, 1) + (s 1 p(:, 2))p(:, 2) 

Evalutating sign2 of this is now a much more interesting affair. If p(:,l) and 
p(:, 2) are orthogonal to one another, i.e., p(:, l) T p(:, 2) = then it is not hard 
to see that both p(:, 1) and p(:, 2) (and their mirrors) will be attractors. In the 
nonorthogonal case we must consider the elements of V = P T P. To get a feel 
for this I encourage you to exercise hop on several 4-by-2 choices. For example, 
with 

/ 1 
1 

-1 

V-i 



p = 



1 \ 

1 
1 

-v 



we find 



V = 



'4 2 N 
2 4 



'p(:,l) T p(:,l) K:>1) T P0,2)\ 
^(:,1) t K: 5 2) p0 5 2) T p(:,2)J " 

and so 1) and p(: : 2) (and their mirrors) are attactors. In fact they are the 
only attractors. A little experimentation reveals that 



[+ + -+] & [+ - - - 

[+- + +], [- + + +] 
[+ - + -], [- + + -] 
[ ], [+--+] 



] 



[+ ] & 

& [+ + + +] 
& [- + -+] 



-> [+ + - -] 

[-- + -] -> [- - + +] 

->[+++ -] 

-> [ +] 
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Project: Face Discrimination via Hopfield 

We will apply Hopfield technology to a variation of last week's task. In par- 
ticular, we will write a function, hopoli, that 

1. Reads and shows the four images, Obama, Clinton, McCain, and Palin and 
converts each to bitstreams of ±1 and displays each as a training pattern as 
in Fig. 14.2. 

2. Constructs and displays the associated Hopfield weight matrix W (as in 
Fig. 14.2) and (iteratively) applies it to very (three fourths of the pixels are 
random) noisy images of our 4 candidates and displays the before and after 
as in Figs. 14.3 and 4. 



Training Pattern 1 Training Pattern 2 The Hopfield Weight Matrix 




Figure 14.2. The black and white candidates and their associated Hopfield 
weight matrix. 



Input Pattern Output Pattern Input Pattern Output Pattern 




Figure 14.3. The men pulled from the noise. 
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Input Pattern Output Pattern Input Pattern Output Pattern 




Figure 14.4. The women pulled from the noise. 

We will use MATLAB's Image Processing Toolbox to read, transform, and show 
these images. In particular, 

a = imread ( ' obama . jpg' ) ; % read the jpg 

abw = im2bw(a); % convert to black and white 

abwc = abw (25 : 435, 60 : end-60 ) ; % crop it 

abwcc = abwc ( 1 : 4 : end, 1 : 4 : end) ; % coarsen it 

imshow (abwcc) % show it 

achieves the first training pattern. This 'final' matrix is now 103-by-82 and is 
of type 'logical' (i.e., T/F or 1/0). To get your bitstream you must reshape 
this into a vector of length 103*82=8446. You must also turn the zeros into 
-ones. I recommend using real. By similar means you can binarize, crop and 
coarsen the other portraits into bitstreams of length 8446. NOTE: This is a 
fairly ambitious bitstream, for it produces a weight matrix W that is 8446-by- 
8446. As MATLAB requires more that 570 megabytes to store this I fear that, 
depending on your hardware, you may receive a memory error. If that is the 
case simply coarsen your image until W fits. 

Your work will be graded as follows: 

6 pts for headers CONTAINING detailed USAGE 
4 pts for further comments in code 
4 pts for indentation 
11 pts for correct W 

11 pts for correct application of W 
14 pts for labeled figures (as above) 
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15. Evolutionary Game Theory 



Please visit the Stanford Encyclopedia of Philosophy for an overview of the 
theory. We will proceed, as usual, by example. In particular, we will play 
Prisoner's Dilemma and the Game of Life in preparation for combination 
of these games in your invesitigation of the Evolution of Cooperation. 

Prisoner's Dilemma 

Two freshmen have been arrested by RUPD for commandeering a golf cart 
and fleeing an officer. They are placed in separate isolation cells. The wily Dean 
makes the following offer to each. 

"You may choose to confess or remain silent. If you confess and your 
accomplice remains silent I will drop all charges against you and charge 
your accomplice with felony theft. Likewise, if your accomplice confesses 
while you remain silent, he will go free while you will be charged with 
felony theft. If you both confess I get two convictions, but I'll see to it 
that you both receive early parole. If you both remain silent, I'll have 
to settle for charging you both with the lesser charge of avoiding arrest. 
If you wish to confess, you must leave a note with the jailer before my 
return." 

The dilemma is that it pays to confess, so long as your accomplice does not also 
confess. 

This game is abstracted as follows. To remain silent is to cooperate with your 
accomplice and so this strategy is denoted by the letter C. To confess is to break 
with, or defect from, your accomplice and so this stategy is denoted by the letter 
D. The payoffs to the row player during a round are 

CD R is the reward for mutual cooperation 

C R S where P is the punishment for mutual defection 
DTP T is the temptation to defect and 

S is for sucker 



and T > R > P > S ensures dilemma 

In this week's assignment you will conduct such games over multiple genera- 
tions and between many many players. In preparation for the incorporation of 
multiple players let us play the 
Game of Life 
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The game is played on an M-by-N board, where the living are ones and the 
dead are zeros, and the rules are brutally simple. 

Any live square with fewer than 2, or more than 3, 

neighbors dies . 
Any dead square with exactly 3 neighbors 

comes to life. 

Please dissect and exercise the associated code. 

Project: Evolution of Cooperation 

We will reproduce the findings of May and Nowak on Evolutionary Games 
and Spatial Chaos. The game takes place on an M-by-N board. Each square is 
occupied by a cooperator (C) or a defector (D). 

In a given round each square plays prisoner's dilemma with each of its imme- 
diate neighbors, including itself. Each square has either 4, 6, or 9 such neighbors. 
The payoffs follows 

C vs . C, each receives 1 

C vs . D, C receives and D receives b 

D vs . D, each receives 

The score of a square is the sum of the neighborly payoffs. 

At the next generation, each square assumes the identity of the winner (high- 
est scorer) of the last neighborly round. 

For example, the score of the board 
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and so if b=1.9 the next generation is 
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D 


D 


c 



We monitor the game's progress by 

painting a square blue when C remains C 
painting a square red when D remains D 
painting a square yellow when C becomes D 
painting a square green when D becomes C 
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You will write the 



function evo (M, N, b, gen) 

that plays this game for gen generations, commencing from a board of co- 
operators with a lone defector at its center (when M and N are odd). With 
M = N = 99 and b = 1.9 you will find boards like 



Generation 10 



Generation 100 



n 




Figure 15.1. The boards at generation 10 and 100. 

after 10 and 100 generations respectively. It is fascinating to watch the fraction 
of cooperators ebb and flow. You will want to quantify this in a graph of the 
form 




100 

Generation 



200 



Figure 15.2. The cooperators take a hit but hold on. 
I recommend that evo delegate to three subfunctions 
function S = score (A, b) 
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function An = advance (S, A) 
function evodisp (A, An) 

where A denotes the state of the previous generation. The state is the matrix 
of identities (C or D, simlilar to live and dead) of each of the players. 

Finally, construct a driver that sets sizes and b in order to generate the plots 
listed below. 

Your work will be graded as follows 

5 pts for headers CONTAINING detailed USAGE 

10 pts for further comments in code 

5 pts for indentation 

20 pts for correct score function 

20 pts for correct advance function 

15 pts for correct evodisp function 

in the following, suppose b=1.9, 

and start from the C board with lone center D 

5 pts for fraction of cooperators plot when M=67, N=67 and gen = 40 
5 pts for fraction of cooperators plot when M=65, N=69 and gen = 100 
5 pts for fraction of cooperators plot when M=69, N=69 and gen = 200 
5 pts for plot of game board at generation 100, when M=199, N=199 
5 pts for plot of game board at generation 200, when M=199, N=199 
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